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Abstract. We study the existence and stability of multisite discrete breathers 
in two prototypical non-square Klein-Gordon lattices, namely a honeycomb and 
a hexagonal one. In the honeycomb case we consider six-site configurations 
and find that for soft potential and positive coupling the out-of-phase breather 
configuration and the charge-two vortex breather are linearly stable, while the 
in-phase and charge-one vortex states are unstable. In the hexagonal lattice, we 
first consider three-site configurations. In the case of soft potential and positive 
coupling, the in-phase configuration is unstable and the charge-one vortex is 
linearly stable. The out-of-phase configuration here is found to always be linearly 
unstable. We then turn to six-site configurations in the hexagonal lattice. The 
stability results in this case are the same as in the six-site configurations in the 
honeycomb lattice. For all configurations in both lattices, the stability results arc 
reversed in the setting of either hard potential or negative coupling. The study is 
complemented by numerical simulations which are in very good agreement with 
the theoretical predictions. Since neither the form of the on-site potential nor 
the sign of the coupling parameter involved have been prescribed, this description 
can accommodate inverse-dispersive systems (e.g., supporting backward waves) 
such as transverse dust-lattice oscillations in dusty plasma (Debye) crystals or 
analogous modes in molecular chains. 
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1. Introduction 

Over the past decade, there has been a considerable increase of interest in the study 
of discrete or quasi-discrete systems in a wide range of areas in physics. Among 
the numerous themes of intense theoretical and experimental interest, one can refer 
to the DNA double-strand dynamics in biophysics [JJ, coupled waveguide arrays 
and photorefractive crystals in nonlinear optics [2j [3j [4], breathing oscillations in 
micromechanical cantilever arrays [5], Bose-Einstein condensates in optical lattices 
in atomic physics [B], granular crystals [7j, and so on. 

These studies have been conducted predominantly in one-dimensional (ID) 
systems or in higher dimensions, but chiefly in the context of square lattices. In 
the latter, a wide variety of novel and interesting phenomena have been revealed, 
both theoretically and experimentally. Pertinent examples include, among others, the 
prediction and observation of dipole 8 , and necklace [9] solitons, discrete vortices 
[10] [TT] , rotary solitons [H [13] , higher-order Bloch modes [14] and gap vortices [15] , 
as well as two-dimensional Bloch oscillations and Landau-Zener tunneling [To] , 

More recently, the dynamics of non-square lattices have become a focal point of 
interest, both in the context of periodic photonic structures [TTJ HH1 EHl [20] , and in that 
of quasi-crystalline [2TJ or completely disordered lattices [22] . While most of the above 
works had a view towards applications based on photorefractive crystals, there exist 
many other applications where such non-square lattices may be relevant. In particular, 
for the hexagonal and honeycomb lattices considered herein, we note that they have 
already been showcased in recent experiments in two-dimensional waveguide arrays 
(e.g., in glass) [23], in optical lattices acting on Bose-Einstein condensates [24], as well 
as in Debye crystals formed in dusty plasmas [25] [26] . Interestingly, in the latter case 
(dusty plasma crystals), the horizontal propagation of transverse (off-plane) vibrations 
gives rise to an inverse-optic dust-lattice mode (a backward wave) [57[ [55J which is 
described by a Klein-Gordon-type equation like the one we shall focus on below, yet 
upon formally considering a negative coupling coefficient ("spring constant"). In this 
case, nonlinearity is provided by the (anharmonic) plasma sheath potential, while 
discreteness is expressed by the (small value of the) ratio of a characteristic coupling 
frequency (related to electrostatic Debye interactions). Remarkably, both ab initio 
theoretical considerations and experimental findings suggest that the sheath on-site 
potential is intrinsically anharmonic and in fact strongly asymmetric, so that existing 
theories involving even (symmetric) polynomial functions for the on-site potential do 
not apply in this case. Further contribution to nonlinearity is furnished by Debye-type 
(screened Coulomb) electrostatic inter-dust-particle interactions, yet the associated 
anharmonicity is of lesser order of magnitude for the transverse mode and can be 
neglected. Details on dusty plasma modelling can be found, e.g., in [29] [30] , while the 
experimental setting is described in 28, 3"T] . 

The above experimental developments have prompted further theoretical work 
towards the goal of understanding the structures that may emerge in such non-square 
lattices and their corresponding stability. An example of this type is the recent work 
of [32] , where this analysis was performed in the framework of the discrete nonlinear 
Schrodinger (DNLS) equation, a prototypical model widely perceived as relevant to 
optical systems. 

The aim of the present work, motivated by models of dusty plasma lattices 
[531 HS]> is t° extend the considerations in [35] to nonlinear Klein-Gordon lattices of the 
hexagonal and honeycomb type. We stress the fact that no assumption will be made on 
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Figure 1. The honeycomb lattice and associated six-site contours 

the (anharmonic) substrate potential or the sign of the coupling parameter involved 
in the description. In particular, inverse-dispersive systems (supporting backward 
propagating waves in the linear limit) such as transverse oscillations in dusty plasma 
crystals are straightforward to accommodate in this description (by reverting, say, the 
sign of the "spring" coupling coefficient to negative) . Within this framework, we offer a 
systematic analysis of the type of solutions that arise in six-site contours in honeycomb 
lattices and three-site as well as six-site contours in hexagonal lattices (for the latter 
lattice, see also earlier work in Refs. 26, 33]). In the six-site contours, we illustrate 
that solutions with higher topological charge S (namely, S = 2) are more dynamically 
robust than lower topological charge ones (namely, S — 1). This is true for soft 
nonlinearities and positive coupling (or hard nonlinearities and negative coupling), 
while the results are reversed if either the nature of the nonlinearity or the sign of the 
coupling are changed. It is interesting to note that similar findings have been reported 
not only in DNLS type chains |32j but also in continuum photorefractive crystals 
|34j (where they have recently been confirmed experimentally 35J). Furthermore, in- 
phase and out-of-phase structures are also examined. The former are found to be 
unstable, while the latter are potentially linearly stable. In the three-site, hexagonal 
configuration, both in- and out-of-phase structures are unstable, while the vortex one 
may be stable. Again, this is true for soft nonlinearities and positive coupling (or hard 
nonlinearities and negative coupling), while all the results, besides the out-of-phase 
case, are reversed if either the nature of the nonlinearity or the sign of the coupling are 
changed. Finally, in the six-site configurations in the hexagonal lattice, the stability 
results are the same as in the case of six-site configurations in the honeycomb lattice. 

Our presentation is structured as follows. In section II, we analyze the six-site 
honeycomb contour, while in section III, we briefly review the findings for the three- 
site hexagonal contour. Section IV contains the comparison with numerical results. 
Finally, section V contains a summary of our findings and a brief discussion of future 
directions. 



2. Existence of multisite breathers in a honeycomb Klein-Gordon lattice 



We consider a honeycomb lattice (see the protoypical example of figure [T]) with on- 
site potential and linear nearest neighbor interaction. This system is described by a 



Existence and stability of multisite breathers in honeycomb and hexagonal lattices 4 



Klein-Gordon Hamiltonian of the form 

h = h q +sh 1 = j2y + v ^ + \ E - ^) 2 ' (!) 

where § is the set of all the oscillators and G is the set of neighboring pairs of 
oscillators. In the anti-continuum (AC) limit 36J (e = 0) we consider the six encircled 
oscillators of figure [1] as moving in periodic orbits with the same period u>i = u>, 
while the rest of the oscillators lie at equilibrium (x i7 pi) — (0,0). This periodic and 
trivially localized motion is continued, for £ / small enough, to provide multi-site 
breathers, if the phase differences between successive oscillators fa are such that they 
correspond to critical points of the effective Hamiltonian H cS 37J. To leading-order 
of approximation, the effective Hamiltonian becomes 

Hf(Ii,<Pi,A) = H (I Z , A) + s{Hi)(Ii,A, fa), 

where Ii,A, fa are canonical variables which are defined by the transformation 



6 

$ = wi, A = J, 



(2) 

fa = w l+ i - urn, Ij = Jj i = 1, ... 5 

j=i+i 

with Ji, Wi being the action-angle variables. The average value of Hi, namely 

(Hi) = — <j) Hidt, is calculated along the unperturbed orbit. Note that, since the 

motion of a single oscillator for e = can be described by Ji = const, and it;, = wt+woi, 
the variables fa can be considered as the phase differences between two successive 
oscillators. Since we consider six oscillators in the AC limit and fa = — Y^i=i &i we 
get only five independent fa. Then, the configurations of the AC limit that will be 
continued to provide the multisite breathers correspond to simple roots of the system 
of equations 



d(H 



— =0, z = l... 5 (3) 



dfa 

provided, of course, the non-resonance of the frequency of the breather ujb with the 
frequency LOb of the system's linear spectrum, namely kujb ^ u) p , and the anharmonicity 
of the single uncoupled oscillator, i.e. dui/dJ ^ 0. This result is in accordance with 
the findings of [38] (see also [33]). 

In order to calculate (Hi) we use the fact that the motion of each oscillator in 
the AC limit can be described by a cosine Fourier series due to the time-reversibility 
[viz. x(-t)=x(t),p(-t) = -p(t)}, 

oo 

x i(t) = X/ ^n(^) cos ( nw i) (4) 

where (J,w) are the action-angle variables for the specific oscillator. Using this fact, 
and the canonical transformation @, (Hi) becomes 

(Hi) = — - X A\ {cosinfa) + cos(nfa) + cos(n0 3 ) + cos(n0 4 ) + cos(n0 5 ) + cos[n(0! + fa + fa + fa\ + fa])}} 
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So, condition ^ becomes 

oo 

^ nA 2 n {sm(n<pi) + sm[n(<f>i + <p 2 + 4>3 + <f>4 + 05)]} = 0, i = 1 ... 5, 

71=1 

which leads to two groups of solutions. The first one is 



0. 



1, 



,5, 



which determines the time- reversible solutions according to the terminology of |36j . 
while the second kind of solutions 



S- 



1,...,5 and Se{l,2} 



determines the non-time- reversible solutions, i.e., vortex breather solutions of 
topological charge S. This kind of solutions can be distinguished from the previous 
one because it possesses a nonzero energy flux |39j . The form and the time evolution 
over a period of the above mentioned motions can be found in 40 . 

The stability of the above mentioned solutions is determined through the 
corresponding Floquet multipliers Xi (see e.g. [41 ). The characteristic exponents 
of the breather are defined as = e aiT . The nonzero characteristic exponents of the 
central oscillators are given as eigenvalues of the stability matrix E = QD 2 H cS where 

SI = ( ^ q ^ is the symplectic structrure matrix and O, I are the 5x5 zero and 

identity matrices respectively. The leading order approximation of E is given by 



/ 



Ex = 



— e 



d 2 (ffi) 



d 2 H t JijH,) 

V oi.i, 



(•5) 



idli I 



xj dlidlj 

If all the eigenvalues of matrix E lie on the imaginary axis, then the corresponding 
breather is linearly stable. If the eigenvalues of E\ are simple to leading-order in s 
and lie on the imaginary axis, then the eigenvalues of E will also lie in the imaginary 
axis, due to continuity reasons. If the eigenvalues of E\ are not simple in the leading- 
order, then the higher-order terms of the approximation could push the eigenvalues 
of E outside the imaginary axis and, thus, lead to complex instability. However, 
this cannot happen if the symplectic signature of E\ is definite [32] and the breather 
remains stable for e small enough. This signature is definite if the quadratic form 
x T E\x has the same sign for every vector x E M 10 . 

Note that, although the above conditions certify linear stability for small values of 
£, for higher values of e a characteristic exponent of the central oscillators can collide 
with the linear spectrum causing instability through a Hamiltonian-Hopf bifurcation, 
as it can be seen in some of our numerical results below (see also |41j). 



2.1. The four configurations under consideration 

We consider four representative configurations of the AC limit which lead to the 
corresponding multisite breathers. For all of the ones we consider, it is true that 
iOi = to <f> J{ = J and (pi = (f>, i = 1...5. More precisely, we consider the in-phase 
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configuration, = 0, the out-of-phase one, cj> 
cf> = 7r/3 and the charge- two vortex breather 
various elements of matrix E\ , 



d 2 H 

BLBL 



= < 



d 2 H, 
' dJ 2 

cPHo 
' dJ 2 




7T, the charge-one vortex breather 
= 27r/3. Using @, we get for the 



duj 

~dJ' 



j = i 

j = i + l 
otherwise 



d 2 (#i) 



f(<P), j ¥= * 



with f(4>) 



\ " „2 42 



2 ^ 

n=l 



n cos(n() 



d 2 (gi) 
d lid I j 



= < 



-2.9i -2g2, i = i + l 
3i, i = * + 2 



, with 



.91 



.92 



1 OO 

-Y 

2 ^ 

n=l 
1 00 



9A» 



2 ^ aj 2 

n— 1 



cos(n0) 



cos(n0) 



OLdd) 




with fc = — n — j A n sm(nq 



n=l 



Below, we examine separately each of the four principal configurations under 
investigation. Up to leading-order of approximation, the eigenvalues of the stability 
matrix E\ correspond to the 0{\fe) approximation of the characteristic exponents, 
and they are given, for all the configurations under consideration, by the expressions: 



C±l,±2 = ± 



dui , 



0~±3±4 



°±5 



It may be added for rigor that we have chosen to keep the sign of the coupling 
strength s arbitrary, bearing in mind that inverse dispersive systems, such as transverse 
dust-lattice vibrations (briefly discussed above), require negative values of e to be 
considered (in contrast with the ordinary Klein-Gordon formulation). Thus, e may be 
either negative or positive throughout this text, unless otherwise stated. 



2.1.1. The in-phase multibreather = 0. We first examine the in-phase (0 = 0) 
configuration. In this case we get dI = 0, which simplifies the calculations of the 
characteristic exponents <7j in equation (|6|). Since for continuous periodic functions 
the size of Fourier coefficients A n is exponentially decreasing as a function of n, f 
converges - in this case to a positive number: 



/ = /(0) 



n=l 



Hence, the nature of the exponents (and the linear stability) hinges on the product ^ 



dJ- 



If it is negative (as is the case for soft nonlinearities and positive couplings, or for hard 
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nonlinearities and negative couplings), the exponents are real and the corresponding 
breather is unstable. If e^jj > we have to check also the corresponding symplectic 
signature to ensure that the exponents will remain in the imaginary axis. It can 
be proved (but the calculations are rather lengthy) that this signature is definite if 
£ TF7 > 0- So, in this case the resulting breather is linearly stable for e small enough 
in order to avoid collisions with the linear spectrum. 

2.1.2. The out-of-phase multibreather, <j> = tt. For 4> = n the stability matrix E\ and 
its eigenvalues (<Tj) are the same as before but in this case it is 



Due to the Hamiltonian nature of the single oscillator it is A\ > AA2 and since the 
absolute value of the terms is exponentially decreasing, / converges to a negative 
number. Following similar arguments as in the previous section we conclude that for 
e§jj < (under the physical conditions discussed in the previous subsection), the 
breather is linearly stable, while it is unstable otherwise. 

2.1.3. The vortex breather of charge S = 1, <j> = 7r/3. In the case of the vortex 

d 2 (H ) 

with charge S — 1 [6 — 7r/3) the terms are not identically zero, but the 

olidfpj 

corresponding characteristic exponents are still, to leading-order of approximation, 
given by equation ([6]). For the polynomial potentials typically used, / converges to a 
positive number: 



Thus, following the arguments of the previous sections for e^j < it turns out that 
the vortex breather is linearly unstable, while for e^j > the breather is linearly 
stable. 

2.1.4- The vortex breather of charge S = 2. The stability matrix E\ and its 
eigenvalues (<7j) are the same as before, but since for = 27r/3, it is found that 



Thus, we conclude that the corresponding breather is linearly stable if e^j < 0, while 
otherwise it is unstable. Note that the symplectic signature arguments mentioned in 
subsection 12.1.11 are still valid. 

3. Existence and stability of 3-site and 6-site breathers in a Hexagonal 
Klein Gordon Lattice 

Let us start by briefly reviewing the theory of Refs. [26l [33], as a preamble to the 
systematic comparison with numerical results in the following section. The Klein- 
Gordon Hamiltonian is of the form of equation fT]), but now each site has six neighbors 
instead of three as before. 



/ = /« = -$> 2 ^(-l)»<0. 





71=1 
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Figure 2. The hexagonal lattice and associated contours. 

We consider an AC limit of three oscillators, as per the contour shown in figure 
[21 In this case, there exist only two independent ^'s since <p 3 — w% — w 2 = — <f>\ — 4>2- 
The persistence conditions for 3-site breathers in these lattices are 

oo 

^2 nA n {sin(n^) + sin [n(<j>i + <fc)]} = i=l,2. 

n=l 

Seeking zeros of the curly bracket in the above expression, we will consider the in-phase 
configuration with <pi = 0, the out-of-phase one with 4>i — tt, and the vortex breather 
with <j>i = =^-. The form and time evolution over a period of these motions can also 
be found in [40) . To leading-order of approximation, the corresponding characteristic 
exponents are given by the following expressions: for the in-phase configuration: 



. •'" jjjl'^h (7) 



for the out-of-phase configuration: 



-n = ±\/-e|j [2/(0) + /(*)], ^=±^-3^/(70, (8) 



and, finally, for the vortex configuration: 



<7±i,± 2 = ±y-3 £ ^/(27r/3), (9) 

where the function f(<j>) is defined as in the previous section. 

In the case of the six-site configurations, the characteristic exponents results are 
the same as in the honeycomb lattice and are given by equations ([5]). This is because 
our theory is a first order one, while the influence of the extra sites in this case would 
be visible in the higher order of the expansion of the exponents. 
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4. Numerical Results 

We perform a set of numerical computations in order to demonstrate the validity of 
our results. Although the theoretical analysis above is completely general, in what 
follows, we consider a particular choice of an on-site potential and breather frequency. 
More specifically, the potential under consideration is a quartic one, of the form, 

V(x) = ^ + b -x* + ^ 

with a = 1, b = —0.27, c = —0.03. For this set of parameters §7 < 0. We will consider 
the orbit with frequency u> = 7.43409 which corresponds to amplitude of oscillation 
x max = 1.949275; thus, J = 1.20306 and §f = -0.224556. For the same orbit we get 
/(0) = 1.423404, /(tt) = -1.279544, / (f ) = 0.638983 and / = -0.710913. 

We point out that a positive value was considered in our numerical investigation 
below. Recalling that the value of e may be either positive or negative (see discussion 
above) , one should therefore keep in mind that the qualitative predictions (on breather 
stability) to follow are directly reversed if either the sign of e or the sign of are 
reversed (to negative/positive, respectively). They obviously remain unchanged if both 
quantities change sign. 

4-1- Honeycomb lattice 

4-1.1. In-phase, 6-site breather. First we consider the in-phase configuration = 
= 0). As mentioned in section !^. 1.11 all the characteristic exponents Oi of the central 
oscillators are real to leading-order of approximation. We calculate the breather 
for increasing values of s and numerically compute the corresponding characteristic 
exponents. In figure [3] we show the theoretically predicted values [cf. equation 
©], depicted by dashed lines, together with the ones calculated by the numerical 
simulation, depicted by solid lines. In this figure, as well as in the following ones, 
shown is only the positive o~i. It is observed that the theoretical and numerical 
branches almost coincide for small values of s, while for larger values of e, where 
the higher-order terms of Oi become significant, the lines start to deviate from each 
other. We note also that, instead of having five branches of numerical o~i, we have 
only three. This happens because the branches of the <7i,2 pair, as well as the 03,4 pair 
coincide (these double pairs will be hereafter denoted by thicker lines), which means 
that the higher-order terms of the approximation also coincide. 

4-1-2. Out-of-phase, 6-site breather. The next configuration is the out-of-phase 
one, i.e. 4>i = cf) = n. As it has already been discussed in section 12.1.21 
this is a linearly stable configuration. Indeed, we have numerically calculated the 
corresponding characteristic exponents for increasing values of the coupling constant 
e and accordingly confirmed that all of them lie on the imaginary axis, as it was 
expected. The calculated value of the exponents, as well as the theoretical 0(y/F) 
prediction, are shown in figure 0] We note again that, for small values of e, the 
theoretical (dashed) with the numerical (solid) curve coincide, and for larger values of 
e, where the higher order terms become significant, they start to diverge. 

4-1.3. Charge S = 1, 6-site vortex breather. We consider now the six-site vortex 
configuration with </>j = = tt/3, i.e., with S — 1. As discussed in section [2.1.31 
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Figure 3. The real part of the characteristic exponents for the in-phase 
honeycomb configuration for increasing values of e. The solid lines represent 
the values obtained numerically, while the dashed ones represent the 0(y/e) 
theoretical prediction. 
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Figure 4. The imaginary part of the characteristic exponents cr; for the out- 
of-phase honeycomb configuration for increasing values of e. The solid lines 
represent the values obtained numerically, while the dashed ones represent the 
0(y/e) theoretical prediction. 

this configuration is unstable. Indeed, the numerical simulation shows that all the 
characteristic exponents of the central oscillators have nonzero real part, as predicted 
by the 0(y/e) estimation of equation ^ [cf. also with the DNLS case of [52] ]. 
These exponents also possess a nonzero imaginary part which is due to higher-order 
contributions of the approximation of <7j. The calculated value of their real part 
with respect to e is shown in the left frame of figure El together with the theoretical 
estimate. Note that the numerical (solid) line and the theoretical (dashed) line diverge 
less than in the two previous cases. In addition in the right frame of figure [5] 
the positive imaginary part of the exponents with the double real part is shown. 
Interestingly, many of these findings (presence of imaginary parts in the next order 
and higher quality of agreement of the real part of the first-order predictions) are also 
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Figure 5. (Left frame) The real part of the characteristic exponents <Tj for 
the charge S = 1 vortex honeycomb configuration for increasing values of e. 
The solid lines represent the values obtained numerically, while the dashed ones 
represent the 0(y/e) theoretical prediction. (Right frame) The imaginary part of 
the exponents with the double real part. 



present in the DNLS equation for these solutions 32J (cf. also the square-lattice case 
in [43]). 

4- 1-4- Charge S — 2, 6-site vortex breather. The configuration with <f>i — <f) — 2tt/3, 
i.e., the S = 2 vortex breather is linearly stable, as it has been already discussed 
in section 12.1.41 Indeed, the numerical simulation confirmed that the characteristic 
exponents of the central oscillators are all imaginary and their value with respect to 
the value of e is shown in figure [HI together with the theoretical prediction. Notice 
that in this case the higher-order splits the double eigenvalue pairs, similarly to the 
case of the DNLS setting [35]. 
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Figure 6. The imaginary part of the characteristic exponents ai for the charge 
5 = 2 vortex breather honeycomb configuration for increasing values of e. The 
solid line represents the values obtained numerically, while the dashed ones 
represent the O(yfe) theoretical prediction. 




4-2. Hexagonal lattice 

4-2.1. In-phase, 3-site breather. First, we consider the in-phase configuration (<pi — 
= 0). The corresponding characteristic exponents are shown in figure [7] The 
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solid line represents the value of the exponent which is acquired from the numerical 
simulation, while the dashed one represents the theoretical prediction acquired from 
equation Q. Note once again in this case the coincidence (also to higher-orders) of 
the double pair of exponents. 
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Figure 7. The real part of the characteristic exponents en for the in-phase three- 
site hexagonal configuration is shown for increasing values of e. The solid line 
represents the values obtained numerically, while the dashed ones represent the 
0(y/e) theoretical prediction. Notice that the double pair does not appear to split 
for the couplings considered. 

4-2.2. Out-of-phase, 3-site breather. The next configuration under consideration is 
the out-of-phase one (fa = <f> = n). This is a linearly unstable configuration, since, 
as seen from equation (jSJ), one of the characteristic exponents will be real while the 
other will be imaginary. This is shown also in panels (a) and (b) of figure [SJ where the 
imaginary and real parts of the exponents are respectively shown. It is observed that, 
for small values of e, one of the exponents is purely imaginary while the other is real. 
When e acquires the value e ~ 0.014 a Hamiltonian-Hopf bifurcation occurs, which 
forms a complex quadruplet of exponents and the corresponding exponent acquires a 
nonzero real part (this collision with a mode of the continuous spectrum leads also to 
the apparent non-smoothness of the line in panel (a)). 

4-2.3. Vortex 3-site breather. The last configuration under consideration is the vortex 
of topological charge S = 1 (fa — 4> = 2ir/S). This configuration is a stable one as it 
can be seen from equation (J9]). In figure [9] the value of the characteristic exponents, 
which are both imaginary, are shown. We see in this case that the two branches 
corresponding to the numerical acquired values of the exponents are distinct, although 
we have only one theoretical prediction (i.e., the double pair splits). This fact implies 
that the higher-order terms contributing to these exponents are different. Notice also 
for £ ~ 0.019, a collision of one of the pairs with the linear spectrum, leading to the 
emergence of an eigenvalue quartet through a Hamiltonian-Hopf bifurcation. Note 
that we do not consider the topological charge S = 2 (fa = <f> = 4tt/3) configuration 
in our study. This happens because the S = 2 case does not produce any physically 
distinct motion, instead it provides the same motion as the 5 = 1 case with the reverse 
rotation direction. 
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Figure 8. The imaginary (a) and the real part (b) of the characteristic exponents 
(T { for the out-of-phase three-site hexagonal configuration for increasing values of 
e. The solid lines represent the values obtained numerically, while the dashed 
ones represent the 0(^/e) theoretical prediction. Notice the Hamiltonian Hopf 
bifurcation at e ~ 0.014, leading to the formation of a quartet. 
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Figure 9. The imaginary part of the characteristic exponents ai for the vortex 
three-site hexagonal configuration for increasing values of e. The solid lines 
represent the values obtained numerically, while the dashed ones represent the 
0(y/e) theoretical prediction. Notice the numerical splitting of the double pair 
and the Hamiltonian-Hopf bifurcation inducing collision for e ~ 0.019. 

4-2.4- In-phase 6-site breather. We consider now 6-site breathers. The first 
configuration we study is the in-phase one (fa — 4> — 0). As it is already predicted this 
configuration is unstable. The corresponding characteristic exponents are shown in 
figure 1101 The solid line represents the value of the exponent which is acquired from 
the numerical simulation, while the dashed one represents the theoretical prediction 
acquired from equation ©. Note once again in this case the coincidence (also to 
higher-orders) of the exponents which have been predicted to be double in the leading 
order of approximation. Although the theoretical estimation of the exponents is the 
same as in the in-phase honeycomb configuration, their actual value is different because 
of the different behaviour of the higher order terms. 

4-2.5. Out-of-phase 6-site breather. The next configuration under consideration is 
the out-of phase one (fa = <j> = 7r), which is linearly stable for small values of e. 
Indeed as it shown in figure [11] all the exponents are purely imaginary until e = 0.011 
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Figure 10. The real part of the characteristic exponents Oi for the in-phase six- 
site hexagonal configuration for increasing values of e. The solid lines represent 
the values obtained numerically, while the dashed ones represent the 0(y/e) 
theoretical prediction. 



where the Hamiltonian Hopf bifurcation occurs. 
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Figure 11. The imaginary part of the characteristic exponents ai for the out-of- 
phase six-site hexagonal configuration for increasing values of e. The solid lines 
represent the values obtained numerically, while the dashed ones represent the 
0(y/e) theoretical prediction. Notice the numerical splitting of the double pair 
and the Hamiltonian-Hopf bifurcation inducing collision for e ~ 0.011. 



4-2.6. S = 1, 6- site vortex breather. We consider now the charge S = 1 vortex 
configuration (fa — <p — 7r/3), which it is anticipated to be linearly unstable, since the 
theoretical prediction (jHJ) implies that all the exponents will be real in 0(y/e). As in 
corresponding case in the honeycomb lattice, the symplectic signature arguments do 
not hold and the exponents, which are predicted to have double real part in leading 
order of approximation, have also a nonzero imaginary part and split, forming two 
complex quartets, for e arbitrary small. This behaviour is shown in figure 1121 In 
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the left frame the real part of the exponents Oi is shown together with the 0{yfe). In 
the right frame the imaginary part of the exponents is shown. Note that, there is no 
theoretical prediction for the imaginary part of the exponents, since this is a higher 
order effect. 




Figure 12. (Left frame) The real part of the characteristic exponents Oi for 
the 5 = 1 vortex six-site hexagonal configuration for increasing values of e. 
The solid lines represent the values obtained numerically, while the dashed ones 
represent the 0(y/e) theoretical prediction. (Right frame) the imaginary part of 
the exponents with the double real part. 



4-2.7. S = 2, 6- site vortex breather. The last configuration under consideration is the 
5 = 2 vortex breather (<fii = (f> = 2tt/3) which is linearly stable. Indeed, as is shown in 
figure [5J the numerical simulation confirmed that the characteristic exponents of the 
central oscillators are all imaginary; their dependence on e is depicted in the figure, 
together with the corresponding theoretical prediction. Notice that, in this case, the 
higher-orders of the characteristic exponents development causes the splitting of those 
which have double imaginary part up to leading order of approximation. 
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Figure 13. The imaginary part of the characteristic exponents <Tj for the 5 = 2 
vortex six-site hexagonal configuration for increasing values of e. The solid lines 
represent the values obtained numerically, while the dashed ones represent the 
0(y/e) theoretical prediction. 
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5. Conclusions 

We have studied the existence and stability of multisite discrete breathers in two 
prototypical examples of non-square two-dimensional lattices. The lattice types under 
consideration were a honeycomb and a hexagonal one. We have considered six- 
site contours in the honeycomb lattice and three-site as well as six-site ones in the 
hexagonal case. We have categorized these solutions in terms of them having zero 
or nonzero energy flux between their bonds, obtaining, respectively, regular discrete 
breather as well as vortex breather structures. 

Our analysis has offered a systematic analysis of the linear stability of the discrete 
excitations supported in the prototypical lattices studied. In the honeycomb lattice 
case we have considered six-site configurations and have shown that, for soft potential 
and positive coupling, the out-of-phase breather configuration and the charge-two 
vortex breather are linearly stable, while the in-phase and charge-one vortex states 
are unstable. In a hexagonal lattice, and in the case of a soft potential and positive 
coupling, the in-phase three-site configuration is unstable and the charge-one vortex is 
linearly stable. The out-of-phase three-site configuration here is always unstable, while 
the stability results for all other configurations in both lattice cases are reversed in 
the setting of either hard potential or negative coupling (but not both). The stability 
results in the six-site hexagonal case coincide with the ones acquired in the honeycomb 
lattice configuration. 

Our study was complemented by numerical computations which have been shown 
to be in very good agreement with theoretical predictions. In order to consider larger 
amplitude multi-breather configurations, we would need a higher-order theory which 
is a natural subject for future investigation. This higher-order treatment could lead us 
also to an improved approximation in the calculation of the characteristic exponents 
considered herein. 

Our results are of relevance in 2D discrete systems characterized by a nonlinear 
substrate potential, namely including molecular chains as well as dusty plasma (Debye) 
lattices. Any particular form/type of on-site potential can fit in our description. 
Finally, we stress that the scope of our formulation includes inverse-dispersive 
(backward wave supporting) systems, such as transverse oscillations in dusty plasma 
crystals. 
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